Probability distribution of the order parameter for the 3D Ising model universality 

class: a high precision Monte Carlo study 
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We study the probability distribution P(M) of the order parameter (average magnetization) M, for 
the finite-size systems at the critical point. The systems under consideration are the 3-dimensional 
Ising model on a simple cubic lattice, and its 3-state generalization known to have remarkably small 
corrections to scaling. Both models are studied in a cubic box with periodic boundary conditions. 
The model with reduced corrections to scaling makes it possible to determine P(M) with unprece- 
dented precision. We also obtain a simple, but remarkably accurate approximate formula describing 
the universal shape of P(M). 
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This work is devoted to the study of the following prob- 
lem. Consider a finite system belonging to the univer- 
sality class of the three-dimensional (3D) Ising model, 
exactly at its critical point. Let the system have a non- 
conserved order parameter, cubic symmetry and periodic 
boundary conditions. For such a finite-size system the 
order parameter M (for the Ising model, the sum of all 
spins, divided by the total number of spins in the system) 
will be a fluctuating quantity, characterized by the prob- 
ability distribution P(M) In the scaling limit (sys- 
tem size going to infinity) this function is universal (up 
to rescaling of M) and can be thus considered a very in- 
teresting and informative characteristic of the given uni- 
versality class. (One should bear in mind that P(M) 
depends on the geometry of the box, and on the bound- 
ary conditions; in this study we always consider a cubic 
box with periodic boundary conditions). For example, 
P(M) contains the information about all momenta (M k ) 
of M, including the universal ratios such as the Binder 
cumulant U = 1 - (l/3)(M 4 )/(Af 2 ) 2 , which has been a 
subject of many Monte Carlo studies Moreover, 
the precise knowledge of P{M) proved to be important 
for locating and characterizing the critical point in Monte 
Carlo studies of various systems, including the liquid-gas 
critical point pd[ |, the critical point in the unified the- 
ory of weak and electromagnetic interactions |l2] ] and in 
quantum chromodynamics | |i3fl . 

The first Monte Carlo computation of P(M) for the 3D 
Ising model in a cubic box with periodic boundary con- 
ditions has been performed in Ref. [Q , where its double- 
peak shape was established. A more accurate determina- 
tion of P(M) has been done in Ref. ]i4| , also by Monte 
Carlo. Results reported for the 3D case in Ref. || ap- 
pear to be incorrect. 

Despite considerable progress in computation of P(M) 
by analytical methods [|l6|-p0|, numerical simulation re- 
mains the main source of information about its proper- 
ties. 



Our aim was to compute P(M) on a qualitatively new 
level of accuracy, in comparison to what has been done 
before ]l4| , and to put the result into form convenient for 
further use. We would like to emphasize the following two 
features of our computation that made this possible. 

1. The computation of Ref. H] used the 3D Ising 
model on a simple cubic lattice of the size 20 3 and 30 3 . 
As we will see, the shape of P(M) obtained with these 
relatively small lattice sizes still differs noticeably from its 
scaling limit, due to nonnegligible corrections to scaling. 
To get over this difficulty, we employed, besides the 3D 
Ising model, the more sophisticated model in the same 
universality class, which was shown to have remarkably 
small corrections to scaling ||. This made it possible to 
determine the scaling limit of P(M) with an accuracy 
far exceeding what would be achievable when one is re- 
stricted to the standard 3D Ising model. 

2. The existing results for P(M) were presented in the 
form of Monte Carlo-generated histograms We 
present a simple 3-parameter formula which is suitable 
for quantitative applications. Its accuracy is about 2 • 
10 -3 of the maximum value of P(M). 

We have performed Monte Carlo simulations of two 
models. The first one is the standard 3D Ising model on 
the simple cubic lattice, defined by the partition function 



y^expj/jy^sj j, Si=±i. 

{Si} 



(1) 



(ij) 



Here (ij) denotes the pairs of nearest neighbours, and the 
sum is over the 2 N possible configurations, where N is 
the total number of spins. We simulate this model at the 
critical point, which we take to be at f3 c — 0.221654, using 
the Swendsen-Wang cluster algorithm [pT[ , an d lattice 
sizes ranging from 12 3 to 58 3 (with periodic boundary 
conditions). 

The second model (with dramatically reduced correc- 
tions to scaling, as was shown in Ref. 0) is the spin-1 
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Blume-Capel model |2^,g3| 
discrete values: —1,0, and 
the partition function 



Here the spins can take 3 
-1. The model is defined by 



Z = ^exp{/3^s i s J - J D^s^} ! Si = -1,0, +1. 

{M (ij) m 

(2) 

The sum thus includes 3 N possible configurations. The 
parameter D is fixed to the special value D = In 2 (as 
explained in Ref. J8|), and we perform the simulations at 
the critical point, which is taken to be (3 C = 0.393422 ||, 
using lattice sizes from 12 3 to 58 3 . The simulations used 
a hybrid method, which alternates one Metropolis sweep 
with 5 or 10 Wolff steps, depending on the system 
size, as described in Ref. ||. 

The probability distribution P(M) is obtained as fol- 
lows. For each configuration generated by the Monte 
Carlo algorithm, we determine the order parameter M = 
j? YliLi Si ' an d increment the population of the corre- 
sponding bin of the histogram by 1. 

We have found that the following ansatz gives a sur- 
prisingly good approximation to P{M): 



r / M 2 \ 2 / M 2 \ i 
P(M) ocexp{-(^-l) (a^ + c)} 



At the same time, the simpler ansatz 



P(M) oc expj- 



-c( 



(3) 



(4) 



is clearly insufficient. This is illustrated in Fig. [j], using 
our highest-statistics data set for the 20 3 lattice. One 
observes that the accuracy of approximation (pT) is ap- 
proximately 20 times higher than that of Eq. (0) , and 
the residual discrepancy of Eq. (|J) is comparable to the 
statistical noise, even with the high statistics used. 

The ansatz (^) was motivated by the observation that 
M 6 plays an important role in the effective potential of 
the models in the 3D Ising universality class, while higher 
powers of the order parameter can usually be neglected 
p5| p?| . That is, the effective potential can in many cases 
well be approximated by a polynomial consisting of M 2 , 
M 4 and M 6 terms. This is exactly what appears in the 
exponent in Eq. (||). 

The approximate nature of the ansatz (^|) manifests 
itself by its failure to correctly reproduce the large-M 
behaviour of the tails of P(M), which is governed by the 
critical index S, 



P(M) oc M (<5_1)/2 cxp{-const ■ M s+1 } 



(5) 



(see Ref. Q; for the discussion of the preexponential 
factor in a more general setting, see Ref. fe6| ). However, 
due to the fact that for the 3D Ising universality class the 
exponent 5 + 1 « 5.8 is close to 6, this does not prevent 
the ansatz ([}]) from accurately describing the main part 
of P{M) (excluding extremely- far-tail region). 



The polynomial in the exponent of Eq. (||) has three 
parameters. Instead of simply parametrizing it by the 
coefficients in front of M 2 , M A and M 6 , we have chosen 
the parametrization so as to separate the scale-invariant 
parameters (a and c) and the scale-dependent parameter 
Mo (which parametrizes the position of the peak of the 
order parameter). The values of a and c in the scaling 
limit are universal and determine the "universal shape" 
of P{M). 

The results of our Monte Carlo simulations are col- 
lected in Tables [| and |l| and shown in Figures [j], ^. For 
the spin-1 model, no deviations from scaling are observed 
on lattices 16 3 and larger, while the simple cubic Ising 
model demonstrates pronounced corrections to scaling, 
which are, even on our largest lattices, much higher than 
both statistical errors of our spin-1 simulations and the 
accuracy of approximation (|jj|). Corrections to scaling 
make it difficult to extract accurate scaling limit values 
of a and c from the simple cubic Ising model data, even 
if statistical errors are reduced by a higher-statistics sim- 
ulation, due to necessity to extrapolate to L — > oo. 

There is no such problem with the spin-1 model, and 
we obtain the universal parameters of Eq. (^) , 



a = 0.158(2), c = 0.776(2). 



(G) 



Here the errors take into account both statistical uncer- 
tainties and the systematic deviations inherent in the ap- 
proximation (^). The latter are estimated from the lower 
right plot in Fig. 

From Eqs. (p|) and @ one can easily obtain any re- 
quired property of P{M). For example, one immediately 
learns that the ratio of the peak value of P(M) to its 
value at M — is 



= 2.173(4). 



(7) 



Summarizing, we have computed, with a higher accuracy 
than previously available, the scaling limit form of the 
probability distribution P{M) of the order parameter M 
of systems with 3D Ising universality, in a cubic box with 
periodic boundary conditions. A convenient description 
of P(M) is given by Eqs. (||) and (||), which deviates 
from the actual P(M) by no more than 2 • 10 -3 times its 
maximum value (Fig. [I], right) . 
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FIG. 1. Probability distribution P(M) of the spin-1 model with reduced corrections to scaling defined by Eq. (Q), and 
its description by approximations given in Eqs. (Q) (left) and ^ (right). Top: P(A-l) obtained by Monte Carlo simulation 
at the critical point (diamonds): 20 3 lattice, /3 = 0.393422, 36 • 10 6 configurations, 1 Metropolis sweep + 5 Wolff steps per 
configuration. The solid line is the best fit with Eq. (|^) (top left) and with Eq. (^) (top right). Bottom: the difference between 
the Monte Carlo data and the fit, corresponding to the plot above it. 
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FIG. 2. Dependence of the scale-invariant parameters a (upper plot) and c (lower plot) of the probability distribution 
P(M), approximated by Eq. (^|), on the lattice size L. The data for the spin-1 model (diamonds) and for the simple cubic Ising 
model (triangles) are taken from Tables | and [n], respectively. The power of the lattice size in the horizontal axis is chosen 
to linearize the leading corrections to scaling, which behave as L - ", where various estimates give uj = 0.80 ± 0.04 (see, e. g., 
Refs. |29l and ffl). 
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TABLE I. The parameters a, c and Md of the probability distribution P(M), approximated by Eq. obtained by the 
Monte Carlo simulation of the spin-1 model defined by Eq. fl2|) at the critical point (/3 = 0.393422). M5W means that a new 
configuration is produced by one Metropolis sweep followed by 5 Wolff steps (see Ref. pf for details) . The last three columns 
are: the scale-invariant (but nonuniversal) quantity MoL d ~ Vh , where yn = 2.4815(15) gj^! X 2 i characterizing the quality of 
fitting the Monte Carlo-generated histogram for P(M) by the ansatz (H), and the number of bins in this histogram. 



Lattice 


Method 


Configs. 


a 


c 


Mo 


AfoL ' 5185 


x 2 


iVbms 


12 3 


M5W 


10 v 


0.1648(15) 


0.7712(19) 


0.30243(16) 


1.0969(6) 


208 


116 


14 :l 


M5W 


3.6 ■ 10 7 


0.1624(9) 


0.7714(9) 


0.27915(9) 


1.0967(4) 


267 


129 


16 3 


M5W 


10 7 


0.1576(18) 


0.7746(14) 


0.26035(11) 


1.0962(5) 


175 


121 


18 :1 


M5W 


10 7 


0.1585(18) 


0.7749(20) 


0.24500(14) 


1.0965(6) 


156 


125 


20 3 


M5W 


9 ■ 10 6 


0.1568(17) 


0.7782(25) 


0.23194(15) 


1.0964(7) 


128 


127 


20 3 


M5W 


3.6 ■ 10 7 


0.1578(8) 


0.7762(12) 


0.23194(7) 


1.0964(3) 


218 


129 


22 3 


M5W 


10 7 


0.1618(16) 


0.7739(21) 


0.22108(11) 


1.0980(6) 


181 


125 


26 3 


M5W 


10 7 


0.1570(26) 


0.7761(27) 


0.20243(14) 


1.0963(8) 


164 


125 


32 3 


M5W 


10 7 


0.1553(19) 


0.7776(24) 


0.18180(11) 


1.0965(7) 


166 


127 


38 3 


M10W 


10 7 


0.1575(16) 


0.7732(25) 


0.16633(10) 


1.0967(7) 


151 


128 


46 3 


M10W 


2 • 10 6 


0.158(5) 


0.781(6) 


0.1506(2) 


1.0964(14) 


136 


123 


58 3 


M10W 


7.2 ■ 10 5 


0.143(7) 


0.776(8) 


0.1331(3) 


1.0927(25) 


135 


119 



TABLE II. Analogous to Table |l| but for the Monte Carlo simulation of the simple cubic Ising model (|l|) at the critical 
point (/J = 0.221654). SW stands for the Swendsen-Wang cluster algorithm pl| . 



Lattice 


Method 


Configs. 


a 


c 


Mo 


MoL 0.5185 


x 2 


-^bins 


12 3 


SW 


7.2 ■ 10 s 


0.268(13) 


0.859(8) 


0.3892(11) 


1.412(4) 


129 


105 


16 3 


SW 


7.2 ■ 10 5 


0.237(6) 


0.845(11) 


0.3360(6) 


1.415(3) 


81.5 


75 


20 3 


SW 


7.2 ■ 10 5 


0.209(9) 


0.839(11) 


0.2984(9) 


1.411(4) 


142 


119 


32 3 


SW 


7.2 • 10 5 


0.200(9) 


0.807(9) 


0.2344(5) 


1.414(3) 


138 


117 


58 3 


SW 


7.2 ■ 10 s 


0.196(14) 


0.807(11) 


0.1733(7) 


1.423(6) 


121 


119 



